Read AquaMaps data
Compiled on Mon Jun 7 15:49:34 2021 by ohara
knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'figs/',
echo = TRUE, message = FALSE, warning = FALSE)
library(raster)
source('https://raw.githubusercontent.com/oharac/src/master/R/common.R')
am_dir <- '~/git-annex/aquamaps_2021'1 Summary
Read in and examine the newest AquaMaps data.
2 Methods
2.1 Read in and examine species summary file
spp_info <- data.table::fread(file.path(am_dir, 'ver10_2019_speciesoccursum_iucn.csv')) %>%
janitor::clean_names() %>%
rename(am_sid = species_id, iucn_sid = iucn_id, comname = f_bname) %>%
mutate(sciname = tolower(paste(genus, species)))Number of unique species (by AM species ID): 33518
Number of unique IUCN species (by IUCN species ID): 10519
2.2 Read in and examine species environmental envelope file
env_info <- data.table::fread(file.path(am_dir, 'ver10_2019_hspen.csv')) %>%
janitor::clean_names() %>%
rename(am_sid = species_id)
names(env_info)## [1] "am_sid" "speccode" "fao_areas"
## [4] "n_most_lat" "s_most_lat" "w_most_long"
## [7] "e_most_long" "depth_yn" "depth_min"
## [10] "depth_pref_min" "depth_pref_max" "depth_max"
## [13] "mean_depth" "pelagic" "temp_yn"
## [16] "temp_min" "temp_pref_min" "temp_pref_max"
## [19] "temp_max" "salinity_yn" "salinity_min"
## [22] "salinity_pref_min" "salinity_pref_max" "salinity_max"
## [25] "prim_prod_yn" "prim_prod_min" "prim_prod_pref_min"
## [28] "prim_prod_pref_max" "prim_prod_max" "ice_con_yn"
## [31] "ice_con_min" "ice_con_pref_min" "ice_con_pref_max"
## [34] "ice_con_max" "oxy_yn" "oxy_min"
## [37] "oxy_pref_min" "oxy_pref_max" "oxy_max"
## [40] "land_dist_yn" "land_dist_min" "land_dist_pref_min"
## [43] "land_dist_pref_max" "land_dist_max" "date_created"
## [46] "layer" "rank" "map_opt"
## [49] "extn_rule_yn"Each parameter (depth, temp, salinity, prim_prod, ice_con, oxy, land_dist) has a column for min, pref_min, pref_max, max (and mean for depth). Additionally, all parameters also have columns for _yn (whether that parameter is relevant for that species). For each parameter, let’s plot distributions of the min, max, and preferred min/max for all spp that depend on that parameter.
params <- c('depth', 'temp', 'salinity', 'prim_prod', 'ice_con', 'oxy', 'land_dist')
prm_or <- paste0(params, collapse = '|')
dist_vals <- c('min', 'pref_min', 'mean', 'pref_max', 'max')
env_sum_df <- env_info %>%
rename(depth_mean = mean_depth) %>%
select(am_sid, contains(params)) %>%
gather(param_full, value, -am_sid) %>%
mutate(param = str_extract(param_full, prm_or),
dist = str_replace(param_full, paste0('(', prm_or, ')_'), '')) %>%
select(-param_full) %>%
spread(dist, value) %>%
filter(yn == 1) %>%
select(-yn) %>%
gather(dist, value, -am_sid, -param) %>%
mutate(dist = factor(dist, levels = dist_vals),
value = ifelse(param == 'depth', log10(value + 1), value),
param = ifelse(param == 'depth', 'log10_depth', param))
write_csv(env_sum_df, file.path(am_dir, 'species_envelope_summary.csv'))
ggplot(env_sum_df) +
geom_density(aes(x = value, color = dist, ..scaled..), fill = NA) +
scale_color_viridis_d() +
facet_wrap(~param, scales = 'free_x') +
theme_minimal()
env_range_df <- env_sum_df %>%
mutate(value = ifelse(param == 'log10_depth', 10^value, value)) %>%
spread(dist, value) %>%
mutate(mid = ifelse(is.na(mean), (pref_min + pref_max)/2, mean),
pref_range = pref_max - pref_min,
max_range = max - min,
pref_range = ifelse(param == 'log10_depth', log10(pref_range), pref_range),
max_range = ifelse(param == 'log10_depth', log10(max_range), max_range)) %>%
arrange(pref_range) %>%
mutate(spp = fct_inorder(am_sid))
ggplot(env_range_df) +
geom_density(aes(x = pref_range, ..scaled..)) +
geom_density(aes(x = max_range, ..scaled..), color = 'darkred') +
facet_wrap(~ param, scales = 'free_x') +
theme_minimal()
2.3 Read in and examine HCAF file
hcaf_info <- data.table::fread(file.path(am_dir, 'hcaf_v7.csv')) %>%
janitor::clean_names() %>%
mutate(across(where(is.numeric), ~ ifelse(.x == -9999, NA, .x)))
names(hcaf_info)## [1] "id" "csquare_code" "loiczid"
## [4] "n_limit" "slimit" "w_limit"
## [7] "e_limit" "center_lat" "center_long"
## [10] "cell_area" "ocean_area" "p_water"
## [13] "clim_zone_code" "fao_area_m" "fao_area_in"
## [16] "country_main" "country_second" "country_third"
## [19] "country_sub_main" "country_sub_second" "country_sub_third"
## [22] "eez" "lme" "lme_border"
## [25] "meow" "ocean_basin" "islands_no"
## [28] "area0_20" "area20_40" "area40_60"
## [31] "area60_80" "area80_100" "area_below100"
## [34] "elevation_min" "elevation_max" "elevation_mean"
## [37] "elevation_sd" "depth_min" "depth_max"
## [40] "depth_mean" "depth_sd" "sst_an_mean"
## [43] "sbt_an_mean" "salinity_mean" "salinity_b_mean"
## [46] "prim_prod_mean" "ice_con_ann" "oxy_mean"
## [49] "oxy_b_mean" "land_dist" "shelf"
## [52] "slope" "abyssal" "tidal_range"
## [55] "coral" "estuary" "seamount"
## [58] "mpa"rast_base <- raster(ext = extent(c(-180, 180, -90, 90)), resolution = 0.5, crs = '+init=epsg:4326') %>%
setValues(1:ncell(.))
var_df <- hcaf_info %>%
select(loiczid, p_water, eez:land_dist,
-starts_with(c('elevation', 'area', 'lme')), -ocean_basin)
vars <- names(var_df)[-1]
for(v in vars) { ### v <- vars[1]
message('plotting map for ', v)
var_rast <- subs(rast_base, var_df, by = 'loiczid', which = v)
plot(var_rast, main = v, col = hcl.colors(n = 100))
}
2.4 Read and examine species cell occurrences
spp_valid <- spp_info %>%
filter(occur_cells >= 10) %>%
select(am_sid, sciname)
# csq_loiczid <- hcaf_info %>%
# select(loiczid, csquare_code) %>%
# distinct()
# spp_cells <- data.table::fread(file.path(am_dir, 'hcaf_species_native.csv')) %>%
# janitor::clean_names() %>%
# left_join(csq_loiczid) %>%
# select(am_sid = species_id, loiczid, prob = probability)
#
# write_csv(spp_cells, file.path(am_dir, 'hcaf_species_native_clean.csv'))
spp_cells <- data.table::fread(file.path(am_dir, 'hcaf_species_native_clean.csv')) %>%
filter(am_sid %in% spp_valid$am_sid)
spp_rich <- spp_cells %>%
group_by(loiczid) %>%
summarize(nspp = n_distinct(am_sid), .groups = 'drop')
spp_rich_rast <- subs(rast_base, spp_rich, by = 'loiczid', which = 'nspp')
plot(log10(spp_rich_rast), col = hcl.colors(20))